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SIMULATIONS  OF  HYPERSONIC  RAREFIED  GAS  FLOWS 

USING  DSMC-MLG 

1  INTRODUCTION 

Under  transitional-flow  conditions,  the  gradients  of  the  macroscopic  variables  become  so 
steep  that  their  scale  length,  L,  is  on  the  order  of  the  average  mean  free  path,  A.  When 
the  Knudsen  number,  Kn  =  X/L,  is  large  (Kn  >  0.03),  the  standard  Navier-Stokes  (NS) 
formulation  cannot  be  used  because  the  constitutive  relations  used  in  the  continuum  formu¬ 
lation  are  not  valid.  For  the  cases  where  Kn  <  0.2,  extensive  efforts  have  been  made  to 
modify  both  the  boundary  conditions  and  the  constitutive  relations  for  the  Navier-Stokes 
equations  (Probstein  and  Pan,  1963;  Oguchi,  1963;  Pan  and  Probstein,  1966;  Kogan,  1969; 
Beskok  and  Karniadakis,  1994)  to  extend  the  range  of  application  of  the  continuum  formu¬ 
lation  to  the  transitional  regime.  However,  the  continuum  approach,  even  with  the  modified 
input  or  formulation,  can  only  be  used  over  a  very  small  range  0.03  <  Kn  <  0.2.  For 
Kn  >  0.2,  it  becomes  necessary  to  use  a  molecular  approach,  such  as  molecular  dynam¬ 
ics  (Allen  and  Tildesley,  1990),  a  Boltzmann  equation  solution  method  (Bhatnagar  et  al., 
1954),  or  a  statistical  molecular  method  such  as  Direct  Simulation  Monte  Carlo  (Bird,  1963). 

The  Direct  Simulation  Monte  Carlo  (DSMC)  method  is  an  approach  that  has  been  used 
widely  and  successfully  to  predict  the  properties  of  transitional-regime  flows.  This  is  a 
numerical  particle-simulation  technique  based  on  kinetic  theory.  The  gas  is  represented  by 
a  large  collection  of  discrete  particles,  which  are  subject  to  intermolecular  collisions  and 
molecule-surface  interactions.  Statistical  accuracy  is  obtained  by  averaging  the  results  over 
many  independent  simulations,  which  may  be  simplified  to  time- averaging  for  steady  flows. 
DSMC  has  been  applied  extensively  to  describe  reentry  flows  (Moss  and  Bird,  1984;  Moss, 
1986;  Bird,  1987;  Moss  et  al.,  1988;  Bird,  1990),  shock  interaction  on  vehicles  at  high 
altitudes  (Carlson  and  Wilmoth,  1992),  and  flows  over  waveriders  (Rault,  1992).  More 
recently,  DSMC  has  been  extended  to  microdynamical  flows,  such  as  microchannels  used  in 
MEMS  devices  (Oh  et  al.,  1995;  Nguyen  et  al.,  1996),  and  to  chemicallyreacting  flows  (Dogra 
et  al.,  1994).  DSMC  has  been  shown  to  be  equivalent  to  solving  the  Boltzmann  equation  for 
problems  involving  monatomic  gases  undergoing  binary  collisions  (Nanbu,  1982). 
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The  DSMC  can  require  very  large  amounts  of  computational  time  and  memory  to  com¬ 
pute  the  properties  of  systems  with  realistic  sizes  and  densities.  The  computational  time 
required  is  directly  proportional  to  N ,  where  N  is  the  number  of  simulated  particles.  Be¬ 
cause  the  statistical  error  reduces  as  1/ y/N. ,  a  large  number  of  samples  must  be  collected  and 
averaged  to  reduce  the  statistical  noise.  The  value  of  N  also  determines  the  grid  resolution. 
For  three-dimensional  flows,  or  flow  with  gradients  that  need  high  local  resolution,  the  use 
of  DSMC  becomes  extremely  limited  by  available  computer  resources.  Recent  attempts  to 
reduce  computational  time  have  involved  both  using  massively  parallel  computers  to  handle 
many  particles  (Furlani  and  Lordi,  1989;  Wilmoth,  1989;  Goldstein  and  Sturtevant,  1989; 
Dagum,  1991)  and  applications  of  new  particle  tracking  methods  (Cybyk  et  al.,  1995). 

In  a  conventional  DSMC  calculation,  the  computational  domain  is  divided  into  a  grid  of 
spatially  fixed  cells.  However,  this  approach  may  result  in  incorrect  prediction  of  collision 
rates  within  sparsely  populated  cells.  In  addition,  the  requirement  of  small  cell  sizes  in 
regions  of  large  macroscopic  gradients  (such  as  shock  structures  and  boundary  layers)  may 
be  violated  at  later  sampling  times.  Both  of  these  problems  are  resolved  by  combining  the 
DSMC  with  the  Monotonic  Lagrangian  Grid  (MLG),  an  algorithm  for  optimizing  particle 
tracking  and  sorting  (Cybyk,  1994).  The  MLG  (Boris,  1986)  is  a  particle  sorting  method  that 
ensures  that  particles  that  are  adjacent  in  physical  space  are  also  adjacent  in  index  space. 
This  means  that,  given  the  location  of  a  particle  in  space  (x,y)  with  array  indices  all  of 
the  nearest  neighbors  of  this  particle  are  offset  by  at  most  a  few  indices  in  any  direction.  This 
is  achieved  by  continually  changing  the  location  of  data  in  computer  memory.  The  method 
is  easily  optimized  on  most  of  types  of  computers,  and  is  extremely  efficient  on  massively 
parallel  computers  (Oh  et  al.,  1996b).  Combining  the  DSMC  and  the  MLG  eliminates  the 
fixed  spatial  grid  of  the  conventional  DSMC  approach,  and  results  in  a  method  that  allows 
automatic  grid  refinement  according  to  the  local  number  density  of  the  gas,  which  in  turn 
produces  higher  accuracy  for  a  given  grid  size  (Cybyk  et  al.,  1995). 

In  this  paper,  the  combined  DSMC-MLG  is  used  on  the  massively  parallel  Connection 
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Machines  to  examine  the  hypersonic,  transitional-regime  flow  through  a  channel  with  a  wedge 
on  the  bottom  surface.  In  the  past,  there  have  been  related  studies  of  hypersonic,  high-Z^n 
flows  on  ramps  and  cones.  For  example,  there  have  been  DSMC  calculations  of  flows  over  a 
two-dimensional  ramp  that  studied  the  pressure  distribution  along  the  ramp  surface  (Moss 
et  al.,  1991),  and  shock/boundary-layer  interactions  as  a  function  of  ramp  angle  (Chpoun 
et  al.,  1993).  A  study  of  flows  over  sharp  cones  (Taylor  et  al.,  1989)  showed  that  the  wake  has 
minimal  effects  on  the  flow  properties  along  the  forebody  of  the  cone,  suggesting  a  similar 
behaviour  in  the  wedge  case.  These  studies  address  parts  of  the  physics  in  the  problem 
analyzed  below. 

Figure  1  is  a  sketch  of  the  typical  flow  that  would  be  expected  for  continuum  conditions. 
Oblique  shocks  start  at  the  leading  edges  of  the  channel  and  the  wedge,  interact  with  each 
other,  and  compress  the  flow  above  the  slanted  part  of  the  wedge  (the  wedge  forebody). 
The  trailing  edge  of  the  wedge  acts  as  a  backward  facing  step,*  so  that  the  flow  expands 
past  the  step  to  extremely  low  density.  The  rarefied  gas  effects,  and  the  inaccuracies  of  the 
NS  formulation,  may  be  most  significant  in  this  region.  Downstream  of  the  wedge,  the  flow 
turns  parallel  to  the  channel  and  forms  a  reattachment  shock.  The  adverse  pressure  gradient 
across  this  shock  causes  the  boundary  layer  to  separate  and  the  flow  reverses  direction. 
Thus  a  recirculating  region  is  created  behind  the  wedge.  These  flow  features-,  effected  by  the 
reflected  and  refracted  shocks,  are  expected  to  change  for  a  rarefied  flow. 

From  a  fluid  dynamic  standpoint,  the  channel-wedge  geometry  produces  a  flow  field  with 
many  important  aerodynamic  features  (such  as  shock  waves,  expansion  fan,  and  bound¬ 
ary  layers)  and  their  complex  interactions,  which  are  not  yet  well  understood  for  rarefied 
conditions.  This  study  shows  the  effects  of  rarefied  conditions  on  the  flow  structures  by  com¬ 
paring  the  computational  results  qualitatively  to  what  would  be  expected  in  a  continuum 
flow  in  Figure  1.  Distributions  of  macroscopic  properties  in  the  flow  field  are  examined.  The 
shear  stress  and  thermal  loading  on  the  solid  surfaces  are  calculated  and  their  relationship 
in  terms  of  the  Reynolds  analogy  is  shown.  The  computed  results  are  then  compared  to 
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existing  theoretical  and  computational  data. 

From  an  algorithmic  standpoint,  the  problem  studied  here  expands  the  application  of 
DSMC-MLG  to  more  complex  geometry  than  the  channel  flows  for  which  it  was  previously 
used.  It  is  thus  possible  to  test  and  develop  aspects  of  the  combined  method,  and  arrive 
at  improved  boundary  conditions  and  more  specific  information  or  the  requirements  for 
convergence  and  the  effects  of  the  quality  of  the  grid  on  the  solution.  In  particular,  we  discuss 
the  problems  of  generating  a  grid  around  the  wedge  and  describe  a  method  for  selecting  the 
best  grid.  Then  we  describe  detailed  resolution  tests  in  which  the  grid  is  successively  refined 
and  the  solutions  are  compared.  We  also  evaluate  what  is  required  to  obtain  statistically 
converged  solution  when  the  time-averaging  process  is  used. 
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2  COMPUTATIONAL  METHOD 


2.1  Direct  Simulation  Monte  Carlo 

The  procedures  for  a  conventional  DSMC  (Bird,  1963)  calculation  are  shown  in  the  flowchart 
in  Figure  2.  The  computational  domain  is  first  divided  into  spatially  fixed  cells  system.  The 
particles  are  then  randomly  distributed  in  the  computational  domain.  These  simulated  par¬ 
ticles,  each  representing  one  or  more  real  gas  molecules,  are  assigned  random  velocities  that 
are  usually  based  on  the  equilibrium  distribution  of  an  undisturbed  gas.  To  begin  the  sim¬ 
ulation,  the  representative  particles  are  moved  for  a  convection  timestep  of  magnitude  A tg. 
This  molecular-motion  process  is  modeled  deterministically.  During  this  step,  interactions 
between  molecules  and  the  solid  boundaries  are  simulated,  and  macroscopic  properties  along 
the  solid  surfaces  can  be  calculated. 

The  next  step  involves  tracking  and  indexing  all  particles.  New  positions  of  the  molecules 
are  sorted.  This  step,  however,  is  a  very  time-consuming  part  of  the  simulation  where 
particles  in  physical  space  must  first  be  indexed  in  computer  memory,  and  their  cell  locations 
must  be  determined  before  collision  pairs  are  selected.  This  step  is  where  the  addition  of  the 
MLG  algorithm  provides  significant  computational  benefits. 

The  molecular-collision  process  is  modeled  probabilistically.  Only  particles  within  a  given 
grid  cell  are  considered  as  possible  collision  partners.  Within  each  cell,  a  representative 
set  of  collisions  is  simulated,  and  collision  pairs  are  selected  randomly.  The  post-collision 
velocities  are  determined  before  particles  are  allowed  to  move  for  the  next  Atg.  This  process 
of  uncoupling  molecular  motions  and  intermolecular  collisions  requires  that  A tg  must  be 
smaller  than  the  mean  collision  time  of  the  unperturbed  gas. 

There  are  serveral  collision  sampling  methods  that  have  been  used  successfully.  Of  those, 
the  time-counter  (TC)  technique  and  the  no-time-counter  (NTC)  technique  (Bird,  1976; 
Bird,  1989)  are  most  commonly  used.  It  has  been  found  that  the  TC  method  allows  the 
acceptance  of  unlikely  collision  pairs  which  results  in  inaccurate  collision  rates.  In  this  study, 
we  use  the  NTC  technique  which  corrects  this  problem. 
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At  selected  time  intervals  At,,  the  macroscopic  properties,  such  as  density,  temperature, 
and  pressure  may  be  sampled.  The  average  of  these  properties  are  calculated  at  the  geometric 
centers  of  the  grid  cells.  When  the  actual-to-simulated  molecule  ratio,  Sm,  is  large,  there 
is  a  high  level  of  statistical  scattering  in  the  results,  and  either  ensemble-averaging  or  time¬ 
averaging  is  needed  to  reduce  the  statistical  fluctuations. 

2.2  Monotonic  Lagrangian  Grid 

The  criteria  ensuring  that  particles  which  are  adjacent  in  physical  space  are  adjacent  in  index 
space  are  the  monotonicity  constraints  (Boris,  1986), 

*(*’»  i)  <  x(i  +  l,j)  for  1  <  *  <  Nx  -  1,  all  j 

y{i,j)  <  y{i,j  +  1)  for  1  <  j  <  Ny  -  1,  all  i,  (1) 

where  Nx  and  Nv  define  the  size  of  the  particle  array  and  x  and  y  are  the  Cartesian  coor¬ 
dinates,  respectively.  A  system  of  particles  that  satisfies  the  monotonicity  conditions  is  said 
to  be  in  MLG  order. 

In  the  types  of  simulations  considered  here,  the  positions  of  the  particles  change  con¬ 
stantly,  leading  to  violation  of  Equation  (1).  The  MLG  order  can  be  restored  by  swapping 
data  stored  in  adjacent  indices  until  the  monotonicity  constraints  of  equation  (1)  are  satis¬ 
fied.  The  data  swapping  is  done  locally.  This  means,  for  instance,  when  particles  move  so 
that  they  are  not  of  MLG  order,  their  MLG  indices  and  the  data  describing  these  particles 
are  exchanged,  or  “swapped,”  until  MLG  order  is  ensured.  Figure  3  shows  the  resorting 
process  for  a  4  x  4  MLG  after  a  convection  time  interval  A tg.  This  resorting  process  results 
in  a.  grid  which  automatically  adapts  to  the  continually  changing  local  number  densities  of 
the  gas. 

The  MLG  subdivides  the  Nx  x  Ny  array  of  simulated  particles  into  nx  x  ny  templates, 
or  blocks  of  nearest-neighbor  particles.  If  nx  =  nv,  the  templates  are  square  in  index  space, 
but  not  necessarily  in  physical  space.  Figure  4  shows  an  example  of  a  nearest-neighbor 
template  that  has  a  population  of  5  x  5  particles.  When  particles  within  an  MLG  template 
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are  allowed  to  interact,  the  template  becomes  analogous  to  the  grid  cell  in  conventional 
DSMC  applications. 

Macroscopic  properties  are  computed  or  “sampled"  at  the  template’s  center-of-mass  and 
depend  on  the  template  area.  In  conventional  DSMC.  the  size  of  the  grid  cells  are  fixed.  In 
DSMC-MLG,  however,  the  continual  change  in  the  size  of  the  templates  requires  continually 
updating  of  the  template  areas.  The  area  of  each  template  is  computed  by  a  very  simple 
and  fast  routine.  Assuming  that  the  Nt  boundary  particles  are  labeled  as  Pi,  P2, ...,  P^b  in  a 
counterclockwise  order,  the  template  area  At  is  given  by 

1  Nt 

A*=o^ ~ xi+i yj)>  (2) 

j=i 

where  xj  and  yj  are  the  coordinates  of  particle  Pj,  and  xsb+ 1  =  Xi,  ysb+ 1  =  y i-  Calculating 
the  template  areas  in  this  way  does  not  account  for  the  area  between  the  template.  To 
remedy  this,  each  template  is  assigned  to  its  own  area  a  fraction  of  the  area  that  is  not 
accounted  for  (Cybyk  et  al.,  1995).  As  a  result,  the  sum  of  all  corrected  template  areas  then 
equals  the  area  of  the  entire  computational  domain. 

2.3  Stochastic  Grid  Restructuring 

For  a  given  system  of  particles,  there  are  many  possible  MLGs  that  satisfy  the  monotonicity 
conditions.  This  is  best  illustrated  by  the  example  in  Figure  5  (Oran  and  Boris,  1987). 
Here,  for  the  same  set  of  16  particles,  three  of  many  possible  grids  are  shown,  all  satisfying 
the  monotonicity  criteria.  As  a  result,  the  quality  of  these  MLGs  can  vary  significantly. 
The  MLG  obtained  by  normal  swapping  of  particle  data  described  above  may  not  be  good 
enough,  while  better  MLGs  often  exist  for  the  same  system  of  particles. 

It  is  possible  to  find  a  high-quality  MLG  for  a  given  spatial  distribution  of  particles  by 
using  Stochastic  Grid  Restructuring  (SGR)  (Sinkovits  et  al.,  1993).  When  this  method  is 
implemented  in  the  DSMC-MLG,  it  is  possible  to  obtain  an  optimal  MLG  for  the  problem 
considered.  The  SGR  consists  of  three  steps: 
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•  The  positions  of  the  particles  are  temporarily  and  randomly  perturbed  by  a  displace¬ 
ment  (6x,6y),  where  —Xd„p  <  Sx  <  xaap  and  — y«fiap  <  Sy  <  ydiap.  During  this 
procedure,  the  unperturbed  positions  are  retained. 

•  The  MLG  swapping  procedures  are  applied  to  the  perturbed  particle  positions  until 
monotonicity  conditions  are  satisfied. 

•  With  the  MLG  for  the  perturbed  positions  as  the  starting  point,  the  MLG  swapping 
procedures  are  applied  to  the  unperturbed  particle  positions  until  monotonicity  condi¬ 
tions  are  satisfied. 

The  quantities  Xdiap  and  ydiap  can  be  chosen  independently,  and  their  magnitudes  greatly 
influence  the  final  MLG.  Several  SGR  iterations  may  be  needed  during  a  simulation  timestep 
to  improve  the  MLG.  Figure  6  illustrates  an  example  of  MLG  improvements  using  SGR.  In 
Figure  6a,  the  MLG  for  a  system  of  6  x  6  particles  without  using  SGR  is  shown.  The  MLGs 
m  Figures  6b  and  c  are  obtained  with  one  and  two  SGR  iterations,  respectively. 

2.4  Parallelization  of  DSMC-MLG 

The  DSMC-MLG  parallelization  method  used  here  is  a  two-level  approach  (Oh  et  al.,  1996b). 
The  first  level  maps  particles  to  an  array  of  processors,  and  the  second  level  maps  templates 
to  an  array  of  processors.  These  arrays  are  of  different  sizes;  there  are  many  more  particles 
than  templates.  This  two-level  approach  allows  efficient  use  of  computer  storage  memory  by 
moving  the  data  at  the  particle-array  level  to  a  much  smaller  array  at  the  template  level, 
thus  reducing  computational  memory  requirements. 

The  flowchart  for  the  parallelized  DSMC-MLG  algorithm  is  shown  in  Figure  7.  The 
shaded  and  striped  blocks  refer  to  the  steps  parallelized  at  the  particle  level  and  template 
level,  respectively.  Particle-level  parallelization  was  performed  for  DSMC-MLG  processes 
that  involve  computations  and  communications  at  the  particle  level  only.  These  were  the 
initialization,  convection,  MLG  resorting,  and  SGR  processes.  Of  these,  the  first  two  steps 
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were  the  simplest,  involving  no  interprocessor  communications.  For  example,  in  the  convec¬ 
tion  step,  all  particles  are  moved  simultaneously  with  their  appropriate  velocities  with  one 
single  instruction.  In  this  step,  some,  but  not  all,  of  the  particles  interact  with  the  bound¬ 
aries.  In  this  case,  the  parallelization  was  implemented  by  selecting  only  those  particles  that 
interact  with  the  boundaries,  performing  the  interactions  in  parallel,  then  updating  their 
velocities  and  positions.  Although  this  process  leads  to  a  load  balancing  problem,  the  effect 
on  the  overall  computing  time  is  minor  (Oh  et  al.,  1996b). 

MLG  resorting  and  Stochastic  Grid  Restructuring  were  also  performed  at  the  particle 
level.  In  DSMC-MLG  applications,  the  particles  move  only  a  small  distance  during  each 
timestep,  and  thus  are  only  slightly  out  of  MLG  order.  The  entire  process  of  resorting  into 
MLG  therefore  requires  only  nearest-neighbor  communication  rather  than  global  communi¬ 
cation.  This  is  ideal,  particularly  on  Connection  Machine  architecture,  for  the  bubble  sort 
algorithm,  which  compares  values  of  nearest-neighbors  of  an  array  on  all  processors  at  the 
same  time. 

The  template  area  (or  cell  volume  for  a  three-dimensional  application),  center  of  mass, 
and  other  properties  are  evaluated  at  the  template  level  for  all  templates  simultaneously. 
The  collision  process  is  also  parallelized  at  the  template  level. 

2.5  Code  Validation 

The  current  code  has  been  validated  by  comparing  the  results  obtained  to  conventional 
DSMC  computations  by  Bird  for  the  Rayleigh  problem  (Cybyk  et  al.,  1995),  the  theoretical 
oblique  shock  angles  for  the  low-Knudsen-number  limit  in  hypersonic  flows  (Oh  et  al.,  1996b; 
Oh  et  al.,  1996a),  and  the  theoretical  Mott-Smith  ratio  of  the  mean  free  path  to  the  shock 
thickness  for  a  monotonic  gas  (Oh  et  al.,  1996a). 

Computational  speed  and  the  DSMC-MLG  have  been  extensively  described  in  the  liter¬ 
ature.  Cybyk  et  al.  (Cybyk  et  al.,  1995)  showed  that  serial  DSMC-MLG  is  indeed  slower 
than  the  conventional  DSMC  approach  in  standard  simplified  test  problems.  However,  there 
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are  significant  advantages  even  in  the  serial  form  for  more  complex  problems,  such  as  those 
that  require  grid  adaptivity  to  resolve  high  gradients  or  moving  bodies.  In  the  DSMC- 
MLG,  there  is,  by  definition,  always  enough  particles  in  each  MLG  cell  to  maintain  optimal 
collision  statistics  (Cybyk  et  al.,  1995).  Oh  et  al.  (Oh  et  al.,  1996b)  showed  very  high  par¬ 
allel  efficiencies  and  speed  ups  of  two  orders  of  magnitude  or  more  compared  to  the  serial 
DSMC-MLG,  and  the  efficiency  increases  dramatically  as  the  number  of  particles  increases. 
DSMC-MLG  computations  are  scalable  and  allow  simulations  of  millions  of  particles  for 
engineering  computation. 

2.6  Inflow-Outflow  Boundary  Condition 

In  the  original  applications  of  the  DSMC-MLG,  the  inflow-outflow  boundary  condition  was 
in  index  space.  The  result  was  that  the  downstream  boundary  could  move  and 
became  distorted  in  physical  space  during  the  simulation  (Cybyk1,  1994).  However,  because 
it  is  most  useful  to  fix  the  length  of  the  system  a  priori ,  a  new  inflow-outflow  boundary 
condition  was  developed  for  this  application  (Nguyen,  1995).  The  procedures  for  maintaining 
this  constant-length  boundary  are  illustrated  in  Figure  8,  which  show  that  the  boundary 
condition  is  applied  in  physical  space  rather  than  in  index  space.  After  the  particles  are 
moved  for  a  timestep  A tg,  those  particles  which  cross  over  a  prescribed  boundary  are  re¬ 
introduced  at  the  inflow  by  resetting  their  velocities  to  freestream  values  and  distributing 
them  at  random  locations  at  the  inflow  boundary.  As  a  result,  the  total  number  of  simulated 
particles  is  kept  constant  for  all  times,  the  total  molecular  mass  in  the  system  is  conserved, 
and  the  length  of  the  system  stays  constant  throught  the  simulation.  Finally,  the  particles 
are  resorted  into  MLG  order.  These  procedures  are  repeated  for  every  convection  timestep. 
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3  GRID,  RESOLUTION,  AND  CONVERGENCE 

STUDIES 


Figure  9  illustrates  the  geometry  and  the  boundary  conditions  for  the  channel-wedge  flow, 
and  Table  1  is  a  summary  of  the  flow  conditions.  The  computational  domain  consists  of 
a  rectangular  channel,  and  a  small  freestream  section  upstream  of  the  leading  edges  of 
the  channel.  A  wedge,  with  half  angle  6 ,  is  placed  on  the  bottom  surface  of  the  channel 
downstream  from  the  leading  edge.  The  channel  is  filled  w'ith  quiescent  rarefied  helium  gas. 
For  the  conditions  in  Table  1,  the  Knudsen  number,  based  on  the  channel  height,  is  0.184, 
which  is  in  the  transitional  flow  regime. 

The  flow  is  initialized  at  Mach  5  everywhere.  The  wall  surfaces  are  kept  at  the  freestream 
temperature  and  are  modeled  as  fully  diffused  boundaries.  The  upper  and  lower  boundaries 
of  the  freestream  section  ahead  of  the  channel  are  specularly  reflecting.  Particles  that  move 
past  the  outflow  boundary  are  reintroduced  at  the  inflow  with  their  properties  reset  to 
freestream  conditions. 

Initial  simulated  results  showed  what  appeared  to  be  a  “blocked”  flow,  with  a  normal 
shock  standing  near  the  inlet  of  the  channel  (Cybyk,  1994).  This  phenomena  in  channel  flow 
is  also  a  common  experimental  problem.  Now  we  initialize  the  calculation  with  molecular 
velocities  that  are  a  small  amount  (3%)  above  those  that  subsequently  flow  into  the  channel. 
This  method  is  very  similar  to  common  experimental  practices  used  to  overcome  choked 
flows  that  often  occur  in  supersonic  and  hypersonic  nozzles. 

3.1  MLG  Grid  Generation  Problem 

The  standard  MLG  sorting  process  does  not  take  into  account  the  presence  of  an  obstacle 
in  the  flow  field.  As  a  result,  in  the  vicinity  of  the  obstacle,  an  MLG  template  could  contain 
particles  from  both  the  front  and  back  sides  of  the  obstacle.  Such  a  template  crosses  the  solid 
boundaries  in  physical  space,  even  though  the  ordering  of  the  particles  satisfies  monotonicity 
conditions.  This  is  illustrated  in  Figure  10,  where  templates  cross  the  wedge  boundaries, 
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and  one  of  the  centers  of  mass  is  inside  the  wedge  and  therefore  outside  of  the  computational 
domain. 

Figure  11  shows  the  evolution  of  the  grid  obtained  by  connecting  the  centers-of-mass  of 
the  templates.  The  presence  of  a  grid  point  inside  the  wedge  corresponds  to  the  presence  of 
a  “saddled”  template,  one  that  saddles  the  wedge.  These  templates  have  extra  area  inside 
the  wedge,  which  in  turn  causes  the  total  area  of  all  the  templates  to  be  larger  than  the 
computational  area.  The  macroscopic  properties  at  the  centers-of-mass  of  these  templates 
are  incorrectly  predicted  as  a  result  of  the  extra  area.  In  addition,  the  particles  in  front 
of  and  behind  the  wedge  are  possible  collision  partners,  which  would  not  occur  in  a  real 
flow.  Therefore,  to  predict  the  flow  properties  accurately,  it  is  necessary  to  make  the  MLG 
conform  to  the  geometry. 

This  problem  is  resolved  by  a  combination  of  two  techniques:  1)  concentration  of  a 
large  number  of  particles  along  the  wedge  surfaces  initially,  and  2)  sorting  locally  in  the 
x-direction  first  before  applying  MLG  sorting.  Figure  12  shows  that  the  resulting  MLGs 
better  conform  to  the  wedge  initially,  and  continue  to  do  so  during  the  simulation.  There 
is  only  one  center-of-mass  inside  the  top-right  corner  of  the  wedge.  However,  the  extra  area 
this  template  adds  is  very  small,  and  the  particles  are  prevented  in  any  event  from  crossing 
the  solid  boundaries.  Any  effect  on  the  final  solution  is  thus  negligible,  and  is  reduced  with 
increasing  grid  resolution. 

This  approach  produced  a  reasonable  grid  for  the  channel-wedge  geometry.  It  was  also 
applied  successfully  on  a  converging-diverging  nozzle  configuration  (Nguyen,  1995).  For  more 
complex  geometries,  however,  it  might  not  give  satisfactory  results,  and  other  techniques 
should  be  used.  In  one  approach,  a  large  number  of  particles  concentrated  on  a  surface 
could  be  used  to  define  the  solid  boundary  and  form  templates.  They  would  not  be  allowed 
to  move  during  the  convection  step,  but  would  be  allowed  to  collide  with  gas  particles  during 
surface  interactions.  Such  extensions  are  currently  being  developed. 
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3.2  Grid  Optimization 

In  Figure  13a,  the  grid  initially  used  for  the  simulation  sometimes  slopes  in  the  direction 
opposite  to  the  flow  structures,  especially  in  the  region  upstream  of  the  wedge.  This  results 
in  a  high  level  of  fluctuations  in  the  solution  (Nguyen,  1995;  Nguyen  et  al.,  1994).  A  better 
MLG  would  result  in  a  smaller  fluctuation  level  and  smoother  results. 

Several  parameters  can  be  used  as  a  measurement  of  the  quality  of  an  MLG.  These  include 
the  average  link  lengths  between  nearest  neighbors  and  the  average  of  the  normalized  dot 
products  of  the  vectors  joining  the  nearest  neighbors  with  the  unit  normals.  The  link  lengths 
are  defined  as 


xunk{i,j)  =  |r(t  +  1,  j)  -  r(i,i)|, 
yiink(i,j)  =  |r (i,j  +  1)  -  r(i,j)|, 


where  r  is  the  vector  joining  the  nearest  neighbors.  The  dot  products  are  defined  as 

_  x(i  +  l,j)-x(i, 

Xlink(l,j) 

.  .x  y(«,i  + 1)  -y(»i 

yunk(i,j) 

For  a  well-structured  MLG,  the  average  link  lengths,  which  indicate  the  average  inter- 
particle  distance,  are  minimized.  The  average  dot  products  are  a  measure  of  the  relative 
direction  between  the  nearest-neighbor  link  and  the  unit  normals,  or  the  orthogonality  of  the 
MLG  grid.  For  example,  if  a  set  of  four  particles  is  arranged  such  that  each  particle  resides 
on  the  lattice  of  a  square,  both  the  x  and  y  dot  products  equal  unity.  A  well-structured 
MLG  has  the  largest  possible  average  dot  products. 

A  series  of  simulations  was  performed  to  examine  the  effects  of  different  MLGs  on  the 
qualities  and  convergence  of  the  computed  flow.  A  nondimensional  SGR  displacement  pa¬ 
rameter,  8 ,  is  defined  as 


8 = 


ydiap 


(3) 
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where  Ar  is  the  averaged  template  size.  Five  different  values  of  S  (0.0,  0.5,  1.0,  2.0,  and  4.0) 
were  considered  for  the  same  grid  resolution  of  1504  templates  (51,144  particles).  Figure  13 
shows  the  MLGs  of  the  steady-state  solutions  obtained  after  2000  timesteps.  There  is  es¬ 
sentially  no  difference  in  the  computing  time  of  the  simulations  with  and  without  the  SGR. 
Once  the  SGR  is  used,  the  slope  of  the  y  grid  reverses  direction  and  is  along  the  direction  of 
flow  structures  (Figures  13b-e).  The  higher  the  value  of  5,  the  more  orthogonal  the  MLG. 

The  average  link  lengths,  x,infc  and  y/infc,  are  shown  in  Figure  14  as  a  function  of  the 
number  of  timesteps  for  all  five  cases.  The  values  of  x/,n*  remain  fairly  constant  during 
the  simulation,  and  are  approximately  the  same  for  all  cases.  The  values  of  ylink  fluctuate 
widely  for  the  case  without  the  SGR.  However,  with  the  SGR,  t //,„*  stays  fairly  constant  with 

time.  When  S  =  2.0,  yunk  is  smallest,  indicating  that  the  average  inter-particle  distances  are 
minimized. 

The  average  dot  products,  xdp  and  ydp,  are  shown  in  Figure  15  as  a  function  of  the  number 
of  timesteps.  The  average  x  dot  products  do  not  vary  much.  However,  the  average  y  dot 
products  fluctuate  widely  without  the  SGR  and  remain  essentially  constant  after  the  first 
few  hundred  timesteps  for  the  other  cases.  Again,  the  MLG  with  5  =  2.0  has  the  highest 
dot  products  and  this  results  in  the  most  orthogonal  grid.  These  results  are  consistent  with 
those  for  the  average  link  lengths. 

Since  without  the  SGR,  and  ydp  fluctuate  widely,  even  after  2000  timesteps,  the 
flow  field  has  high  level  of  statistical  scattering,  and  is  unfit  for  time-averaging.  An  optimal 
MLG  appears  to  be  the  case  where  S  =  2.0.  However,  it  is  unclear  which  grid  gives  the  best 
solution  for  this  problem  since  grid  orthogonality  is  not  required  in  these  calculations.  When 
S  =  0.5,  the  MLG  is  highly  adapted  to  the  local  number  density  of  the  flow.  As  S  increases, 
the  orthogonality  of  the  MLG  increases,  but  it  appears  to  conform  less  to  the  flow  structure. 
An  answer  was  found  in  the  density  contours  in  Figure  16  that  compares  the  solutions  of 
the  highest  quality  MLG  (£  =  2.0)  to  the  MLG  most  closely  resembling  the  flow  (<S  =  0.5). 
The  flow  fields  were  obtained  by  averaging  500  timesteps  after  the  steady  state  was  reached 
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after  2000  timesteps.  Even  though  the  results  appear  somewhat  similar,  a  closer  comparison 
show  that  for  8  =  0.5,  the  solution  is  better  resolved,  especially  in  the  regions  of  the  oblique 
shocks  and  their  interaction  at  the  channel  centerline.  Figure  17  shows  that  the  density 
and  temperature  along  the  centerline  of  the  channel  are  indeed  very  similar,  except  that  the 
8  =  0.5  case  shows  somewhat  higher  values  of  density  in  the  shock-interaction  region.  The 
analysis  of  the  flow  field  solution,  presented  in  a  later  section,  is  for  the  8  =  0.5  case. 

3.3  Grid  Resolution  Study 

A  series  of  simulations  of  the  channel-wedge  problem  was  performed  to  study  the  sensitivity 
of  the  results  to  the  number  of  MLG  templates.  Table  2  gives  the  properties  of  four  different 
grids,  with  the  number  of  templates  ranging  from  900  to  1800  and  the  number  of  simulated 
particles  ranging  from  32,400  to  64,800.  In  all  of  these  cases,  the  SGR  parameter  8  is  2.0. 
The  analysis  could  also  be  performed  for  8  =  0.5,  the  most  adaptive  grid  case.  However,  as 
seen  in  the  previous  section,  results  from  the  two  cases  are  essentially  the  same,  and  therefore 
no  significant  effect  on  the  analysis  should  be  expected. 

The  grids  for  the  four  cases  are  shown  in  Figure  18  and  the  corresponding  density  contours 
are  shown  in  Figure  19.  The  flow  features  obtained  using  the  coarsest  grid  are  not  nearly  as 
well  resolved  as  those  using  1504  and  1800  templates,  which  are  very  comparable.  Figures  20 
show  the  density  distributions  along  the  bottom  and  top  surfaces  of  the  channel.  The  results 
are  comparable  for  all  cases  except  the  coarsest  resolution.  In  particular,  the  differences 
between  the  1504-  and  1800-template  cases  are  small,  indicating  grid-independence  of  the 
solution.  Because  a  simulation  on  the  1504-template  grid  requires  35%  less  computing  time 
than  that  on  the  1800-template  grid,  the  1504-template  grid  is  used  for  further  studies. 

3.4  Time- Averaging  and  Convergence 

Due  to  the  high  ratio  of  actual  to  simulated  molecules,  there  are  large  statistical  fluctuations 
which  must  be  greatly  reduced  to  obtain  a  meaningful  solution.  This  is  be  accomplished  by 
time-averaging  the  solution  once  the  steady  state  is  established.  Using  the  1504-template 
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grid  with  S  =  0.5,  the  flow  reaches  steady  state  after  about  1000  timesteps.  The  simulation 
was  then  continued  for  a  few  hundred  additional  timesteps  before  time-averaging  was  started. 

Figure  21  shows  the  density  contours  averaged  over  40,  100,  300,  and  500  timesteps 
beyond  steady  state.  This  is  equivalent  to  averaging  the  solution  over  40,  100,  300,  and  500 
independent  ensembles,  respectively.  As  the  number  of  averaged  timesteps  (or  ensembles) 
increases,  the  statistical  fluctuations  decrease.  The  results  for  the  300  and  500-timestep 
averagings  are  practically  the  same,  indicating  convergence. 

Since  the  location  and  size  of  the  MLG  templates  may  change  with  time,  it  is  interesting 
to  consider  whether  the  grid  moves,  and  how  much,  during  the  time-averaging  process.  The 
calculations  showed  that  at  approximately  700  timesteps  into  the  pre-averaged  solutions, 
the  grid  no  longer  changed  significantly.  This  is  in  agreement  with  the  results  described  in 
Figures  14  and  15,  which  show  the  properties  of  the  MLG  —  the  averaged  link  lengths  and 
dot  products  —  remain  constant  with  time  after  the  few  hundred  initial  timesteps.  The 
MLGs  taken  at  various  sampling  times,  shown  in  Figure  22,  are  basically  identical. 

Figure  23  compares  the  density  along  the  top  and  bottom  surfaces  of  the  channel  for 
time-averaging  stages  of  40,  300,  and  500  timesteps.  The  distributions  are  essentially  the 
same,  with  the  40-timestep  case  showing  more  fluctuations. 
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4  ANALYSIS  OF  THE  FLOW  FIELD 


4.1  Flow  Features  and  Surface  Properties 

Figure  24  shows  the  converged  solution,  represented  by  density,  temperature,  pressure,  and 
Mach  number  contours,  for  the  1504-template  grid  with  S  =  0.5.  The  flow'  field  is  similar  to 
that  in  Figure  1,  only  now  all  of  the  flow  structures  are  very  diffused.  However,  it  is  still 
possible  to  identify  many  features. 

The  density  contours  in  Figure  24a  show  the  boundary-layer  growth  along  the  solid 
surfaces.  The  thickness  of  the  boundary  layers  increases  as  the  flow  moves  downstream. 
Near  the  outflow  boundary,  the  flow  expands  as  it  leaves  the  channel,  causing  the  boundary- 
layer  thickness  to  decrease.  The  temperature  contours  in  Figure  24b  show  the  temperature 
of  the  gas  in  the  core  of  the  channel  is  much  higher  than  that  near  the  walls,  which  are  held 
at  a  constant  temperature  of  273  K.  As  a  supersonic  flow  expends  through  a  rarefraction 
fan,  the  density,  temperature,  and  pressure  of  the  gas  decrease.  This  effect  is  seen  in  the 
region  immediately  downstream  of  the  trailing  edge  of  the  wedge  in  Figures  24a  —  c. 

The  pressure  contours  in  Figure  24c  show  two  oblique  shocks  emanating  from  the  upper 
and  lower  leading  edges  of  the  channel  and  intersecting  at  about  10  cm  downstream  of  the 
leading  edges.  The  oblique  shock  angle  is  ~  22°.  This  was  estimated  by  measuring  the 
angle  of  the  pressure  contours  at  the  leading  edge  of  the  channel.  The  pressure  contours  also 
show  that  there  is  a  third  oblique  shock  at  an  angle  of  ~  40°  emanating  from  the  leading 
edge  of  the  wedge.  All  of  these  shock  structures  interact  near  the  centerline  of  the  channel 
above  the  wedge  forebody  in  a  way  similar  to  that  shown  in  Figure  1.  In  the  upper-half  of 
the  channel,  a  refracted  shock  interacts  with  the  boundary  layer  near  the  upper  wail  and 
reflects  downstream.  In  the  lower-half  of  the  channel,  a  weaker  refracted  shock  merges  with 
the  rarefraction  region  that  was  caused  by  the  supersonic  expansion  of  the  flow  around  the 
sharp  corner  of  the  trailing  edge  of  the  wedge.  The  flow  then  eventually  turns  parallel  to  the 
wall  to  form  a  very  diffuse  viscous  layer  downstream.  In  a  continuum  flow  as  in  Figure  1,  this 
layer  would  be  a  reattachment  shock.  The  adverse  pressure  gradient  (Figure  24c)  across  this 
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reattachment  region  causes  the  boundary  layer  to  separate,  and  a  reverse  flow  is  developed 
just  behind  the  trailing  edge  of  the  wedge.  The  pressure  reaches  a  local  maximum  at  the 
bottom  wall  near  x  =  33  cm  as  the  flow  crosses  the  reattachment  layer.  The  Mach  contours  in 
Figure  24d  show  a  supersonic  core  flow  between  the  subsonic  boundary  layers  and  expansion 
near  the  outflow  boundary. 

Figures  25  shows  the  density,  pressure,  temperature,  and  magnitude  of  velocity  along  the 
top  and  bottom  solid  surfaces  of  the  channel.  Along  the  bottom  surface,  at  about  x  =  5 
cm,  the  pressure  rises  sharply  across  the  leading-edge  shock  and  increases  significantly  along 
the  wedge  surface  due  to  the  presence  of  the  second  oblique  shock.  Just  ahead  of  the  wedge 
trailing  edge,  the  pressure  drops  across  the  rarefraction  region  as  the  flow  expands,  then 
increases  again  as  the  boundary  layer  develops.  Along  the  top  surface,  the  pressure  increases 
across  the  leading-edge  shock  in  the  same  manner  as  along  the  bottom  surface.  The  pressure 
then  rises  markedly  at  about  x  =  15  cm  due  to  the  interactions  between  the  boundary  layer 
and  the  refracted  shock. 

Figure  25c  shows  the  temperature  of  the  gas  adjacent  to  the  walls.  From  this  we  see  that 
the  temperature  jump  is  about  four  times  the  wall  temperature  near  the  leading  edge  and 
decreases  rapidly  further  downstream.  The  magnitude  of  velocity  in  Figure  25d  is  what  is 
called  the  velocity  slip  along  the  walls.  At  the  entrance  to  the  channel,  the  slip  velocity  is 
high  because  the  density  is  low.  As  the  boundary  layer  develops,  downstream,  the  density 
increases  and  the  slip  velocity  drops  sharply  to  near  zero.  A  small  increase  in  the  slip  velocity 
near  the  exit  plane,  x  =  40  cm,  indicates  that  the  flow  expands  as  it  leaves  the  channel. 

The  development  of  the  viscous  boundary  layers  is  shown  in  Figure  26  by  the  velocity 
profiles  at  various  locations  along  the  x-axis.  At  the  leading  edge  of  the  channel,  the  bound¬ 
ary  layers  begin  to  form,  and  they  grow  as  the  flow  moves  downstream.  The  velocities  are 
negative  at  the  wall  for  20  cm  <  x  <  22  cm,  indicating  a  reverse  flow. 

The  development  of  the  thermal  boundary  layers  is  shown  in  Figure  27  by  the  temperature 
profiles  at  various  locations  along  the  x-axis.  Compared  to  the  velocity  profiles  in  Figure  26, 
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the  temperature  profiles  are  not  as  symmetric,  reflecting  different  shock  strengths  due  to  the 
presence  of  the  wedge.  The  temperature  profiles  also  show  the  temperature  jump  near  the 
walls. 


4.2  Comparisons  of  Flow  Properties 


For  an  inviscid  flow,  with  no  boundary  layers  and  infinitely  thin  walls,  the  leading-edge 
shocks  degenerate  to  Mach  waves.  For  Moo  =  5  flow,  the  angle  of  such  a  Mach  wave  would 
be  11.54°,  which  is  much  smaller  than  the  22°  of  the  shock  angle  calculated  in  this  study, 
indicating  the  displacement  effect  of  the  boundary  layer.  The  boundary  layer  grows  rapidly 
in  the  leading  edge  region  and  acts  as  an  effective  body  that  deflects  the  incoming  streamlines 
upward.  As  a  result,  the  shock  at  the  leading  edge  is  nearly  normal  to  the  surface  of  the 
channel.  Across  this  shock  at  the  surface,  P2/P1  is  approximately  5.  This  is  much  smaller 
than  the  pressure  rise  of  29  across  a  continuum  normal  shock,  and  is  a  consequence  of  the 
low  density  conditions. 

When  a  high-velocity  flow  is  slowed  by  the  presence  of  a  surface  boundary,  the  viscous 
dissipation  in  the  boundary  layer  transforms  the  high  kinetic  energy  of  the  molecules  into 
internal  energy  of  the  gas.  This  causes  a  large  increase  in  the  gas  temperature,  which  in 
turn  causes  the  viscosity  coefficient  to  increase.  The  net  effect  is  that  the  boundary  layer 
becomes  thicker  and  grows  more  rapidly  than  it  would  for  lower-speed  flows  at  the  same 
Reynolds  number  (Anderson,  1989).  This  boundary  layer  causes  a  major  displacement  effect 
on  the  outer  inviscid  flow.  In  turn,  the  changes  in  the  inviscid  flow  affect  the  properties  of 
the  boundary  layer.  This  phenomenon  is  called  viscous  interaction. 

In  continuum  flow,  the  viscous  interaction  near  the  leading  edge  may  be  described  by  the 
similarity  parameter  defined  as  (Anderson,  1989) 
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The  quantity  C  is  the  Chapman-Rubesin  parameter  and  is  defined  as 


where  the  subscript  w  indicates  wall  properties.  In  the  continuum  flow  regime,  \  dictates 
the  induced  pressure  along  the  surface  of  a  flat  plate.  For  a  cold-wall  flat  plate,  the  induced 
pressure  for  a  weak  interaction  is  (Hayes  and  Probstein,  1959) 

~  =  1  +  0.078*.  (6) 

In  low-density  flows,  the  pressure  interaction  near  the  leading  edge  should  be  less  pro¬ 
nounced  than  that  in  the  continuum  flow  due  to  slip  effects.  Figure  28  compares  the  pressure 
near  the  leading  edge  of  the  top  surface  of  the  channel  with  that  calculated  by  Equation  (6). 
The  pressure  in  the  low-density  case  is  consistently  smaller  than  that  in  the  continuum  ap¬ 
proximation.  Also,  the  pressure  drop  in  the  vicinity  of  the  leading  edge  (the  region  indicated 
by  the  smaller  values  of  1  /{Ax)  in  the  abscissa)  is  due  to  the  fact  that  the  leading-edge  shock 
diffuses  out  over  the  leading  edge  of  the  channel. 

Figure  29  compares  the  results  from  the  present  study  to  those  obtained  from  a  DSMC 
calculation  (Chpoun  et  al.,  1993)  for  a  near-continuum,  =  4  flow  over  a  20°  compression 
corner.  Figure  29a  shows  the  surface  pressure  distribution  as  a  function  of  the  nondimen- 
sionalized  distance  x/ L j,  where  Lj  is  5  cm,  the  length  of  the  flat  section  upstream  of  the 
leading  edge  of  the  wedge.  Due  to  the  higher  freestream  Mach  number,  the  pressure  in  the 
present  study  is  consistently  higher  than  that  from  Chpoun  et  al. 

Figure  29b  compares  the  velocity  distribution  along  the  ramp  surface.  The  more  rarefied 
condition  ( Kn  =  0.184)  in  the  present  work  resulted  in  a  more  gradual  slip  in  the  velocity 
adjacent  to  the  surface.  In  contrast,  the  near-continuum  condition  (Kn  =  0.0047)  in  Chpoun 
et  al.  caused  a  much  sharper  decrease  near  the  leading  edge  of  the  flat  section. 

4.3  Skin  Friction  and  Heat  Transfer 

The  skin  friction  on  the  walls  of  the  channel  is  calculated  by  averaging  the  parallel  momentum 
transfer  of  the  molecules  impinging  on  the  surfaces.  Let  the  parallel  momentum  of  an  incident 
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molecule  be  Pino  and  that  of  the  same  molecule  after  reflection  be  Pre/i,  such  that 


Pnc  =  mMline’ 

Pre/l  =  mMlr ejr  (") 

where  m  is  the  mass  of  each  molecule,  and  Vjj  is  the  component  of  velocity  parallel  to  the 
surface.  The  total  change  in  the  parallel  momentum  of  the  molecule,  AP||,  is  given  by 

AP||  =  Pinc  -  Preft.  (8) 


The  skin  friction,  tw,  on  a  small  surface  area,  AA,  over  an  incremental  time,  A tg,  generated 
by  the  impinging  molecules  is 

_  sM 

TW  ~  A  A  \  A  > 


AAAt0 


(9) 


where  the  summation  sign  indicates  the  total  contribution  from  all  impinging  molecules. 
During  a  DSMC  calculation,  the  skin  friction  on  the  surfaces  of  the  solid  wails  may  be 
sampled  at  any  time  using  Equation  (9).  The  heat  flux,  qw,  at  the  walls  may  be  calculated 
in  a  similar  manner  by  replacing  the  parallel  momentum  with  the  energy  of  the  impinging 
molecule. 

The  results  are  presented  in  terms  of  the  dimensionless  skin  friction  coefficient,  c/,  and 
the  Stanton  number,  st,  which  are  given  by  (Anderson,  1989) 


Cf  ~  1  .  t/2  ’ 


I 


(10) 


and 


Qw 

St  =  -  A„)  ’  (11> 

where  the  subscript  e  refers  to  conditions  at  the  edge  of  the  boundary  layer,  hw  is  the 
enthalpy  at  the  wall,  and  haw  is  the  enthalpy  if  the  wall  is  assumed  to  be  adiabatic.  The 
adiabatic  wall  enthalpy  is  evaluated  in  terms  of  the  recovery  factor,  r,  by  (Anderson,  1989) 

V 2 

hatu  =  he  -f-  r  —  .  (12) 
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The  recovery  factor  is  taken  equal  to  y/Pr ,  where  Pr  is  the  Prandtl  number  (Anderson, 
1989).  For  helium,  Pr  =  0.7  (White,  1974).  In  this  study,  the  freestream  conditions  are 
assumed  at  the  edge  of  the  boundary  layers. 

The  skin  friction  and  heat  transfer  to  the  top  and  bottom  surfaces  of  the  channel  are 
shown  in  Figure  30.  Along  the  bottom  surface,  the  value  of  cj  is  highest  near  the  leading 
edge  where  the  high-velocity  molecules  are  slowed  down  substantially  by  the  solid  boundaries. 
The  skin  friction  coefficient  decreases  downstream,  then  increases  along  the  wedge  surface, 
reaching  a  local  maximum  near  the  end  of  the  wedge.  Behind  the  trailing  edge  of  the  wedge, 
cj  becomes  negative  in  the  recirculation  region.  Futher  downstream,  c/  stays  fairly  constant, 
with  a  slight  increase  near  the  exit  plane  where  the  flow  expands. 

Along  the  top  surface,  cj  peaks  near  the  entrance  plane  and  decreases  downstream.  A 
local  maximum  of  cj  occurs  near  x  =  25  cm  as  a  result  of  the  interaction  between  the  refracted 
shock  and  the  boundary  layer  adjacent  to  the  top  wall.  The  distributions  of  st  are  similar  to 
those  of  cf.  The  Stanton  number  is  positive  everywhere,  including  the  recirculating  region, 
indicating  the  transfer  of  heat  from  the  flow  to  the  walls.  Figure  30  represents  an  average  of 
1000  timesteps.  These  figures  indicate  that  the  skin  friction  and  heat  transfer  require  more 
timestep  averaging  than  the  state  properties  to  obtain  the  same  level  of  fluctuations. 

The  results  along  the  top  wall  agree  well  with  those  obtained  by  a  DSMC  calculation 
of  a  hypersonic  flow  along  a  flat  plate  (Moss  et  al.,  1991).  In  their  results,  both  the  skin 
friction  and  heat  transfer  peak  near  the  leading  edge  and  decrease  downstream.  The  peak  cj 
in  Moss  et  al.  (1991)  is  approximately  0.07  compared  to  0.056  for  the  top  wall  in  this  study. 
The  peak  st  in  Moss  et  al.  (1991)  is  0.02,  whereas  it  is  about  0.03  for  the  top  wall.  The 
quantitative  differences  are  due  to  the  different  freestream  and  wall  temperature  conditions. 

In  incompressible  flow  theory,  the  Stanton  number  and  the  skin  friction  coefficient  along 
a  laminar  flat  plate  are  related  by  the  Reynolds  analogy  (Anderson,  1989), 

1  _2 

st  =  c/—  Pr  3.  (13) 

Even  though,  the  flow  is  hypersonic  in  this  study,  it  is  interesting  to  see  how  much  the  results 


22 


in  Figure  30  deviate  from  the  above  relation.  Figure  31  compares  two  distributions  of  the 
Stanton  number:  one  from  the  DSMC  simulations  (Figure  30),  and  the  other  calculated  by 
substituting  the  computed  cj  into  the  Reynolds  analogy  of  Equation  (13).  In  Figure  31a, 
the  two  distributions  of  st  along  the  top  wall  compare  very  well,  except  for  the  region  near 
the  interaction  between  the  refracted  shock  and  the  boundary  layer.  In  Figure  31b,  the 
two  distributions  are  about  the  same  along  the  flat  section  of  the  bottom  wall  ahead  of 
the  wedge.  Downstream  of  this  section,  however,  the  Reynolds  analogy  values  deviate  from 
the  DSMC  values  due  to  the  presence  of  the  wedge.  The  results  in  these  figures  show  that 
the  Reynolds  analogy  might  still  be  applicable  for  high- A/,  high- if n  flow  along  flat  plates 
with  the  freestream  and  wall  temperature  conditions  similar  to  those  in  this  study.  Other 
freestream  and  wall  conditions  must  be  investigated  before  the  Reynolds  analogy  can  be 
generalized. 
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5  CONCLUSION 


Simulations  of  a  rarefied,  high-A’n  flow  of  helium  gas  through  a  channel  with  a  wedge  have 
been  performed  using  the  combined  DSMC-MLG  algorithm.  These  simulations  were  used 
both  to  describe  the  flow  features  and  their  interactions  under  rarefied  conditions,  and  to 
develop  and  test  a  new  approach  to  the  DSMC  algorithm.  Because  this  work  had  a  fluid 
dynamic  as  well  as  algorithmic  objective,  an  extensive  series  of  computations  was  performed 
for  a  single  physical  system. 

The  resulting  computations  present  a  highly  resolved  picture  of  a  flow  field  that  was 
compared  qualitatively  to  what  would  be  expected  in  a  continuum  flow.  Where  possible,  the 
results  were  compared  quantitatively  to  theoretical  analysis  and  computations.  Since  the 
flow  is  so  rarefied,  the  flow  features,  such  as  oblique  shocks,  rarefraction  fan,  and  boundary 
layers,  and  their  interactions  are  diffuse.  The  pressure  jump  across  the  leading-edge  shock 
at  the  surface  is  considerably  less  than  that  for  a  continuum  normal  shock.  Comparisons 
to  hypersonic  weak  viscous  interaction  theory  for  continuum  flow  show  that  the  pressure 
interaction  near  the  leading  edge  is  less  pronounced  under  rarefied  conditions.  Downstream 
of  the  trailing  edge  of  the  wedge,  the  reattachment  shock  that  is  typical  in  a  continuum  flow 
degenerates  into  a  very  diffuse  viscous  layer. 

The  effects  of  low-density  conditions  are  also  shown  in  the  slip  in  velocity  and  temperature 
of  the  gas  adjacent  to  the  solid  surfaces.  The  slip  velocity  is  highest  at  the  entrance  of 
the  channel  where  the  density  is  low,  and  decreases  more  gradually  than  that  for  a  near- 
continuum  conditions  in  the  work  of  Chpoun  et  al.  (1993).  Similarly,  the  temperature  jump 
peaks  near  the  leading  edges  reaching  approximately  four  times  the  wall  temperature,  and 
then  decreases  further  downstream  as  the  boundary  layers  develop. 

The  skin  friction  and  heat  transfer  calculations  show  that  peak  shear  force  and  thermal 
loading  occur  near  the  leading  edges  of  the  channel.  The  effects  of  the  wedge  on  the  flow 
field  are  shown  in  local  maxima  of  skin  friction  coefficient  and  Stanton  number  distributions 
on  the  bottom  wall.  In  addition,  viscous  interactions  between  the  refracted  shock  and  the 
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boundary  layer  near  the  top  wall  result  in  an  increase  in  skin  friction  and  heat  transfer.  These 
calculations  compare  fairly  well  with  DSMC  calculations  on  a  flat  plate  by  Moss  et  al.  (1991). 
The  Stanton  number  distribution  along  the  top  wall  compare  well  with  that  calculated  by 
the  Reynolds  analogy  for  incompressible  boundary-layer  flow,  suggesting  that  the  Reynolds 
analogy  might  be  applicable  for  high-speed,  rarefied  flow  along  flat  plates  under  conditions 
similar  to  those  in  this  study.  These  results,  along  with  the  slip  velocity  and  temperature 
jump  data,  provide  a  benchmark  for  testing  phenomenological  and  theoretical  models  that 
attempt  to  extend  the  applicability  of  the  NS  equations  to  other  flow  regimes. 

Another  important  aspect  in  this  study  is  the  development  and  implementation  of  a  new 
inflow-outflow  boundary  condition.  Previous  applications  of  DSMC-MLG  used  an  inflow- 
outflow  boundary  condition  that  causes  the  downstream  boundary  to  move  in  physical  space 
during  the  simulation.  The  new  boundary  condition  effectively  keeps  the  length  of  the 
the  channel-wedge  geometry  constant  throughout  the  simulation  while  conserving  the  total 
molecular  mass  in  the  system. 

The  MLG  provides  the  DSMC  with  automatic  grid  adaptation  according  to  local  number 
density  and  gives  correct  collision  rates  everywhere  in  the  flow  field.  The  DSMC-MLG 
combination,  when  implemented  on  a  massively  parallel  computer,  is  extremely  efficient  and 
allows  simulations  of  very  large  numbers  of  particles.  For  the  specific  problem  of  the  flow 
through  a  channel- wedge,  a  typical  calculations  of  65,000  particles  took  about  20  minutes  or 
less  to  converge  on  the  256  node  CM-5.  There  were  a  number  of  practical  questions  about 
using  the  DSMC-MLG  that  have  now  been  resolved.  This  work  has  provided  substantial 
insight  into  issues  of  the  grid  generation  in  the  presence  of  an  obstacle,  the  resolution  required 
in  terms  of  the  number  of  simulated  particles,  the  effects  of  grid  structure  on  the  solution 
and  how  to  find  an  optimal  grid,  and  what  is  required  in  terms  of  averaging  the  solution  to 
eliminate  the  effects  of  statistical  fluctuations. 
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Table  1:  Flow  conditions  for  the  channel- wedge  simulation. 


Definition  Flow  Condition  Value 

Parameter 


length  of  channel 
height  of  channel 
wedge  half-angle 
freestream  density 
freestream  temperature 
freestream  Mach  number 
wall  temperature 
freestream  mean  free  path 
Knudsen  number 
number  of  particles/template 


%max 

35  cm 

Umax 

7.5  cm 

e 

10° 

Pug 

2.0e-9  g/cm: 

Too 

273  K 

Moo 

5 

Tw 

273  K 

XUg 

1.38  cm 

Kfl  =  XUg  / Umax 

0.184 

Nt 

6x6 

Table  2:  Computational  parameters  for  the  grid  sensitivity  study. 


Number  of 
Templates 

Grid  Size 

Molecules 
per  Template 

Total  Number  of 
Simulated  Molecules 

75  x  12 

6x6 

32400 

6x6 

43200 

94  x  16 

6x6 

54144 

100  x  18 

6x6 

64800 

30 


leading-edge 

shocks 


reattachment 


shock 


Figure  1:  A  typical  continuum  flow  field  in  the  channel-wedge  geometry. 
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START 


AVERAGE  RUNS  FOR  STEADY  FLOW  or 
AVERAGE  SAMPLES  TAKEN  AFTER 
ESTABLISHMENT  OF  STEADY  FLOW 


/  PRINT  FINAL  RESULTS  / 

^  1  — J 

dsTOPj> 

Figure  2:  Flowchart  for  a  conventional  Direct  Simulation  Monte  Carlo  computation. 
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a)  timet 


b)  time  t  +  Atg, 
before  resorting 


c)  time  t  +  Atg, 
after  resorting 


Figure  3:  Example  of  an  MLG  resorting  process. 
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Figure  4:  Example  of  a  5x5  template  of  nearest  neighbors. 
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Figure  5:  An  example  of  three  grids  that  are  consistent  with  the  same  MLG  monotonicity 
constraints  [Oran  and  Boris,  1987] 
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SET  CONSTANT 


CONSTRUCT  MLG 
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<?>  H«  ’  fc/ftaS#  h'ft'i’V-  Vi  t-)  ;i  *)  h  ; ; 


Figure  7:  A  flowchart  of  the  parallelized  DSMC-MLG  algorithm  [Oh  et  al.,  1995] 
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a) 


c) 


Figure  8:  The  procedures  used  in  the  constant-length  boundary  condition,  a)  The  particles 
are  move  for  a  A tg,  b)  those  that  move  past  x  =  xmax  are  re-introduced  at  the  inflow,  c) 
these  particles  are  assigned  freestream  conditions. 
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Figure  9:  Schematic  of  the  channel-wedge  flow  numerical  experiment. 
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Figure  10:  Schematic  showing  how  templates  may  cross  solid  boundaries. 


Figure  11:  MLGs  for  the  channel-wedge  simulation.  Grids  do  not  conform  to  the  geometry, 
a)  Initial  grid,  b)  after  50  timesteps,  c)  after  100  timesteps. 
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Figure  12:  MLGs  for  the  channel- wedge  simulation.  Grids  conform  to  the  geometry,  a) 
Initial  grid,  b)  after  50  timesteps,  c)  after  100  timesteps. 
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gure  13:  MLGs  of  the  channel- wedge  flow  problem  using  SGR  with  five  different  S  values. 
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Figure  15:  Average  dot  products  for  t] 
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Figure  16:  Comparsion  of  density  contours  of  the  channel-wedge  flow  field  for  different 
displacement  parameters. 
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Figure  17:  Comparsion  of  computed  density  and  temperature  distributions  along  the  cen¬ 
terline  of  the  channel. 
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Figure  18:  MLGs  for  the  grid  resolution  study. 
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Figure  19:  Density  contours  for  the  four  grid  resolutions. 
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Figure  20.  Comparison  of  density  distributions  for  tbe  four  grid  resolutions; 
b)  bottom  surface  of  the  channel. 
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Figure  21:  Flow  field  density  contours  for  a)  40,  b)  100,  c)  300,  and  d)  500  time-averaging 
steps. 
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Figure  22:  Comparison  of  1504-template  MLGs  for  a)  40,  b)  100,  and  c)  500  time-averaging 
steps 
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Figure  25:  Macroscopic  properties  along  the  channel-wedge  solid  surfaces:  a)  density,  b) 
pressure,  c)  temperature,  d)  magnitude  of  velocity. 
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Figure  27:  Temperature  profiles  at  various  x-locations  along  the  channel. 
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1/(Ax) 

Figure  28:  Comparison  of  calculated  normalized  pressure  distribution  near  the  top  leading 
edge  of  the  channel  to  that  of  weak  viscous  interaction  theory.  A  =  0.5(7  —  1)0.664(1  + 
2.6TW/T0). 
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Figure  29:  Comparisons  of  surface  properties  with  data  from  Chpoun  et  al.,  (1993).  a) 
Pressure,  b)  Velocity. 
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Figure  31:  Comparisons  of  DSMC-calculated  and  Reynolds  analogy  Stanton  number  distri¬ 
butions:  a)  top  surface,  b)  bottom  surface. 
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